!*************************************************************
! INTERPOLATION PROCEDURES
!*************************************************************
 module interp
    implicit none
    integer, parameter:: prec=8
    real(prec), parameter:: one=1.0_prec, zero=0.0_prec
    
  private

  !******* followings can be called from outside *******
  public:: lip
  public:: blip
  public:: spline
  public:: splint  
  public:: splie2
  public:: splin2
  public:: splie2_lip_spl
  public:: splin2_lip_spl
  public:: lininterp3
    contains
! **************************************************

  ! subroutine: LININTERP
        SUBROUTINE LININTERP(fx,x,x1,x2,fx1,fx2)
        IMPLICIT NONE
        real(prec) :: fx, x,x1,x2,fx1,fx2
        fx=((x-x2)/(x1-x2))*fx1+((x-x1)/(x2-x1))*fx2
        RETURN
        END SUBROUTINE LININTERP

!***************************
! LINEAR INTERPOLATIONS
!***************************
!*********************************************************************
!24apr13
  function lip_old(n,xvec,zmat,x)
  ! x value to evaluate; xvec(n) variable vector; zmat - function to evaluate
	integer :: n
	real(prec) :: xvec(1:n), zmat(1:n), x
    real(prec) lip_old
	integer :: i, j
	real(prec) ::  t, u
! Check for invalid input
	if ((x.lt.xvec(1)) .or. (x.gt.xvec(n))) then
		write(*,*)  'x outside of grid! Program stops,BLIP'
		stop
	endif
! Locate the interval that surrounds (x)
	i=count(xvec.le.x)
! Determine return value
	if ((i.eq.n)) then
		lip_old=zmat(n)
	else
		t=(x-xvec(i))/(xvec(i+1)-xvec(i))
		lip_old=(one-t)*zmat(i)+t*zmat(i+1)
	endif
  end function lip_old
  
!*********************************************************************
  function lip(xvec,zmat,x)
  ! x value to evaluate; xvec(n) variable vector; zmat(n) - function to evaluate
	integer :: n
	real(prec),dimension(:) :: xvec, zmat
    real(prec) :: x
    real(prec) lip
	integer :: i, j
	real(prec) ::  t, u
    n=size(xvec)
! Check for invalid input
	if ((x.lt.xvec(1)) .or. (x.gt.xvec(n))) then
		write(*,*)  'x outside of grid! Program stops, LIP'
		!stop
		write(*,*)  'x ', x,  xvec(1), xvec(n)
	endif
! Locate the interval that surrounds (x)
	i=count(xvec.le.x)
! Determine return value
	if ((i.eq.n)) then
		lip=zmat(n)
	else
		t=(x-xvec(i))/(xvec(i+1)-xvec(i))
		lip=(one-t)*zmat(i)+t*zmat(i+1)
	endif
  end function lip


! BILINEAR INTERPOLATIONS
!*********************************************************************
!12apr13
  function blip_old(n,m,xvec,yvec,zmat,x,y)
	integer :: n, m
	real(prec),dimension(:) :: xvec(1:n), yvec(1:m), zmat(1:n,1:m)
    real(prec) :: x, y, blip_old
	integer :: i, j
	real(prec) ::  t, u
! Check for invalid input
	if ((x.lt.xvec(1)) .or. (x.gt.xvec(n))) then
		write(*,*)  'x outside of grid! Program stops,BLIP'
		stop
	endif
	if ((y.lt.yvec(1)) .or. (y.gt.yvec(m))) then
		write(*,*)  'y outside of grid! Program stops,BLIP'
		stop
	endif
! Locate the square that surrounds (x,y)
	i=count(xvec.le.x)
	j=count(yvec.le.y)
! Determine return value
	if ((i.eq.n) .and. (j.eq.m)) then
		blip_old=zmat(n,m)
	else if ((i.eq.n) .and. (j.lt.m)) then
		u=(y-yvec(j))/(yvec(j+1)-yvec(j))
		blip_old=(one-u)*zmat(n,j)+u*zmat(n,j+1)
	else if ((i.lt.n) .and. (j.eq.m)) then
		t=(x-xvec(i))/(xvec(i+1)-xvec(i))
		blip_old=t*zmat(i+1,m)+(one-t)*zmat(i,m)
	else
		t=(x-xvec(i))/(xvec(i+1)-xvec(i))
		u=(y-yvec(j))/(yvec(j+1)-yvec(j))
		blip_old=(one-t)*(one-u)*zmat(i,j)+t*(one-u)*zmat(i+1,j)+t*u*zmat(i+1,j+1)+(one-t)*u*zmat(i,j+1)
	endif
  end function blip_old
!*********************************************************************
  function blip(xvec,yvec,zmat,x,y)
	integer :: n, m
	real(prec),dimension(:) :: xvec(:), yvec(:), zmat(:,:)
    real(prec) :: x, y
    real(prec) ::blip
	integer :: i, j
	real(prec) ::  t, u
    n=size(xvec)
    m=size(yvec)
! Check for invalid input
	if ((x.lt.xvec(1)) .or. (x.gt.xvec(n))) then
		write(*,*)  'x outside of grid! Program stops,BLIP'
		write(*,*)  'x ', x,  xvec(1), xvec(n)
!		stop
	endif
	if ((y.lt.yvec(1)) .or. (y.gt.yvec(m))) then
		write(*,*)  'y outside of grid! Program stops,BLIP'
		write(*,*)  'y ', y,  yvec(1), yvec(n)
!		stop
	endif
! Locate the square that surrounds (x,y)
	i=count(xvec.le.x)
	j=count(yvec.le.y)
! Determine return value
	if ((i.eq.n) .and. (j.eq.m)) then
		blip=zmat(n,m)
	else if ((i.eq.n) .and. (j.lt.m)) then
		u=(y-yvec(j))/(yvec(j+1)-yvec(j))
		blip=(one-u)*zmat(n,j)+u*zmat(n,j+1)
	else if ((i.lt.n) .and. (j.eq.m)) then
		t=(x-xvec(i))/(xvec(i+1)-xvec(i))
		blip=t*zmat(i+1,m)+(one-t)*zmat(i,m)
	else
		t=(x-xvec(i))/(xvec(i+1)-xvec(i))
		u=(y-yvec(j))/(yvec(j+1)-yvec(j))
		blip=(one-t)*(one-u)*zmat(i,j)+t*(one-u)*zmat(i+1,j)+t*u*zmat(i+1,j+1)+(one-t)*u*zmat(i,j+1)
	endif
  end function blip

!*********************************************************************
!	3 dimensional linear interpolation
!	using continuous convex combination of neighboring grid points 
! Note: adapted on July 2014 from Tony Smith (search 'chged')

 function lininterp3(xgrid, ygrid, zgrid, ZZ, x, y, z) !chged

	implicit none

	integer :: nx, ny, nz, xndx, yndx, zndx
	real(prec) :: t, u, v, xyz
	real(prec), intent(in) :: x, y, z, xgrid(:), ygrid(:), zgrid(:), ZZ(:,:,:)
    real(prec) :: lininterp3
    
	nx = size(ZZ,1)
	ny = size(ZZ,2)
	nz = size(ZZ,3)

	xndx = locate(xgrid, x) !chged
	yndx = locate(ygrid, y)
	zndx = locate(zgrid, z)

	if (xndx == 0 .or. xndx == nx) then
		print*, "the first argment is out of range in lininterp3"
	else if (yndx == 0 .or. yndx == ny) then
		print*, "the second argment is out of range in lininterp3"
	else if (zndx == 0 .or. zndx == nz) then 
		print*, "the third argment is out of range in lininterp3"
	end if

	t = (x - xgrid(xndx))/(xgrid(xndx+1) - xgrid(xndx))
	u = (y - ygrid(yndx))/(ygrid(yndx+1) - ygrid(yndx))
	v = (z - zgrid(zndx))/(zgrid(zndx+1) - zgrid(zndx))

	xyz = t*u*v*ZZ(xndx+1, yndx+1, zndx+1)				&
		+ (1-t)*u*v*ZZ(xndx, yndx+1, zndx+1)			&
		+ t*(1-u)*v*ZZ(xndx+1, yndx, zndx+1)			&
		+ (1-t)*(1-u)*v*ZZ(xndx, yndx, zndx+1)			&
		+ t*u*(1-v)*ZZ(xndx+1, yndx+1, zndx)			&
		+ (1-t)*u*(1-v)*ZZ(xndx, yndx+1, zndx)			&
		+ t*(1-u)*(1-v)*ZZ(xndx+1, yndx, zndx)			&
		+ (1-t)*(1-u)*(1-v)*ZZ(xndx, yndx, zndx)

	lininterp3 = xyz

end function



! ***************************************
! CUBIC SPLINES MODULES
! ***************************************
!11dec12
  subroutine spline(x,y,yp1,ypn,y2)
    implicit none
    real(prec), dimension(:) :: x, y
    real(prec) :: yp1, ypn
    real(prec), dimension(:) :: y2
     !Given tabulated arrays x and y of length n, 
     ! and 1st derivatives yp1 ypn, it returns array 
     ! y2 of 2nd derivatives
     
    integer :: n
    real(prec), dimension(size(x)) :: a, b, c, r
    n=size(x)
    
    c(1:n-1)=x(2:n)-x(1:n-1)
    r(1:n-1)=6.0_prec*((y(2:n)-y(1:n-1))/c(1:n-1))
    r(2:n-1)=r(2:n-1)-r(1:n-2)
    a(2:n-1)=c(1:n-2)
    b(2:n-1)=2.0_prec*(c(2:n-1)+a(2:n-1))
    b(1)=1.0
    b(n)=1.0
    if(yp1.gt.0.99E+30_prec) then
        r(1)=0.0
        c(1)=0.0
    else
        r(1)=(3.0_prec/(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
        c(1)=0.5
    end if
    if(ypn.gt.0.99E+30_prec) then
        r(n)=0.0
        a(n)=0.0
    else
        r(n)=(-3.0_prec/(x(n)-x(n-1)))*((y(n)-y(n-1))/(x(n)-x(n-1))-ypn)
        a(n)=0.5       
    end if
    call tridag(a(2:n),b(1:n),c(1:n-1),r(1:n),y2(1:n))    
  end subroutine spline    
  !************************************************    
  subroutine tridag(a,b,c,r,u)
  implicit none
  real(prec), dimension(:) :: a,b,c,r,u
  real(prec), dimension(size(b)) :: gam
  integer :: n, j
  real(prec) :: bet
  
  n=size(b)
  bet=b(1)
  if(bet.eq.0.0) then
    write(*,*)' Error in tridag st1 !! '
    pause
  end if
  
  u(1)=r(1)/bet
  do j=2, n
   gam(j)=c(j-1)/bet
   bet=b(j)-a(j-1)*gam(j)
   if(bet.eq.0.0) then
    write(*,*)' Error in tridag st1 !! '
    pause
   end if
   u(j)=(r(j)-a(j-1)*u(j-1))/bet
  end do
  do j=n-1, 1, -1
   u(j)=u(j)-gam(j+1)*u(j+1)
  end do   
  end subroutine tridag    
 !********************************************
 function splint(xa,ya,y2a,x)
    implicit none
    real(prec), dimension(:) :: xa, ya, y2a
    real(prec) :: x, splint
     !xa, ya - tabuated arrays
     !y2a - aray output from spline
     !x - value at which splint gives interpolated value 
    integer :: khi, klo, n
    real(prec) :: a,b,h
    n=size(xa)
    klo=max(min(locate(xa,x),n-1),1)
    
    khi=klo+1
    h=xa(khi)-xa(klo)
    if(h==0.0) then
     write(*,*)' ad xa input in splint st1 !! '
     pause    
    end if
    a=(xa(khi)-x)/h
    b=(x-xa(klo))/h
    splint=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**2)/6.0_prec
 end function splint
 !********************************************
 function locate(xx,x)
    implicit none
    real(prec), dimension(:) :: xx
    real(prec) :: x
    integer :: locate
    integer :: n, jl, ju, jm
    logical :: ascnd
    n=size(xx)
    ascnd= (xx(n) >= xx(1))
        jl=0
        ju=n+1
        do
         if(ju-jl<=1) exit
         jm=(ju+jl)/2
         if(ascnd.eqv.(x >= xx(jm))) then
            jl=jm
         else
            ju=jm
         end if
        end do
        if(x==xx(1)) then
         locate=1
        else if (x==xx(n)) then
         locate=n-1
        else
         locate = jl
        end if
 end function locate

! 2 DIM SPLINES
 !********************************************
 subroutine splie2(x1a,x2a,ya,y2a)
    implicit none
    real(prec), dimension(:) :: x1a, x2a
    real(prec), dimension(:,:) ::  ya, y2a
    integer :: j,m,ndum
     !Given m x n function ya, ... (from Press et al)
     ! x1a firt variable; x2a second varaible 
     ! ya tabulated function
     ! y2a 2nd derivatives of second variable
    m=size(x1a)
    ndum=size(x2a)     
    
    do j=1,m
     call spline(x2a,ya(j,:),1.0e30_prec,1.0e30_prec,y2a(j,:))
    end do
 end subroutine splie2

 !********************************************
 function splin2(x1a,x2a,ya,y2a,x1,x2)
    implicit none
    real(prec), dimension(:) :: x1a, x2a
    real(prec), dimension(:,:) ::     ya, y2a
    real(prec) :: x1, x2, splin2
    integer :: j,m,ndum
     !Given m x n function ya, ... (from Press et al)
     ! x1a firt variable; x2a second varaible 
     ! ya tabulated function
     ! y2a 2nd derivatives of first variable from splie2
     ! x1, x2 - values where to evaluate by bicubic interpolation
    real(prec), dimension(size(x1a)) :: yytmp, y2tmp2
    m=size(x1a)
    ndum=size(x2a)     
    
    do j=1,m
     yytmp(j)=splint(x2a,ya(j,:),y2a(j,:),x2)
    end do
     call spline(x1a,yytmp,1.0e30_prec,1.0e30_prec,y2tmp2)
     splin2=splint(x1a,yytmp,y2tmp2,x1)
 end function splin2


 !******END CUBIC SPLINES **************************

!**********************
! 2 DIM SPLINES+LINEAR
!**********************
!********************************************
 subroutine splie2_spl_lip(x1a,x2a,ya,y1a)
    implicit none
    real(prec), dimension(:) :: x1a, x2a
    real(prec), dimension(:,:) ::  ya, y1a
    integer :: j,m,ndum
     !Given m x n function ya, ... (from Press et al)
     ! x1a firt variable; x2a second varaible 
     ! ya tabulated function
     ! y1a 2nd derivatives of first variable
    m=size(x2a)
    ndum=size(x1a)     
    
    do j=1,m
     call spline(x1a,ya(:,j),1.0e30_prec,1.0e30_prec,y1a(:,j))
    end do
 end subroutine splie2_spl_lip

 !********************************************
 function splin2_spl_lip(x1a,x2a,ya,y1a,x1,x2)
    implicit none
    real(prec), dimension(:) :: x1a, x2a
    real(prec), dimension(:,:) ::     ya, y1a
    real(prec) :: x1, x2, splin2_spl_lip
    integer :: j,m,ndum
     !Given m x n function ya, ... (from Press et al)
     ! x1a firt variable; x2a second varaible 
     ! ya tabulated function
     ! y1a 2nd derivatives of 2nd variable from splie2
     ! x1, x2 - values where to evaluate by cubic+linear interpolation
    real(prec), dimension(size(x2a)) :: yytmp 
    m=size(x2a)
    ndum=size(x1a)     
    
    do j=1,m
     yytmp(j)=splint(x1a,ya(:,j),y1a(:,j),x1)
    end do
    
    splin2_spl_lip=lip(x2a,yytmp,x2)

 end function splin2_spl_lip

!********************************************
 subroutine splie2_lip_spl(x1a,x2a,ya,y2a)
    implicit none
    real(prec), dimension(:) :: x1a, x2a
    real(prec), dimension(:,:) ::  ya, y2a
    integer :: j,m,ndum
     !Given m x n function ya, ... (from Press et al)
     ! x1a firt variable; x2a second varaible 
     ! ya tabulated function
     ! y2a 2nd derivatives of 2nd first variable
    m=size(x1a)
    ndum=size(x2a)     
    
    do j=1,m
     call spline(x2a,ya(j,:),1.0e30_prec,1.0e30_prec,y2a(j,:))
    end do
 end subroutine splie2_lip_spl




!********************************************
 function splin2_lip_spl(x1a,x2a,ya,y2a,x1,x2)
    implicit none
    real(prec), dimension(:) :: x1a, x2a
    real(prec), dimension(:,:) ::     ya, y2a
    real(prec) :: x1, x2, splin2_lip_spl
    integer :: j,m,ndum
     !Given m x n function ya, ... (from Press et al)
     ! x1a firt variable; x2a second varaible 
     ! ya tabulated function
     ! y1a 2nd derivatives of 2nd variable from splie2
     ! x1, x2 - values where to evaluate by linear+spline interpolation
    real(prec), dimension(size(x1a)) :: yytmp 
    m=size(x1a)
    ndum=size(x2a)     
    
    do j=1,m
     yytmp(j)=splint(x2a,ya(j,:),y2a(j,:),x2)
    end do
    
    splin2_lip_spl=lip(x1a,yytmp,x1)

 end function splin2_lip_spl
!********************************************
 end module interp
